This document demonstrates how to create various visualizations of
the EQ5D instrument using ggplot2 in
R.
It also serves as a user manual for an accompanying clean R
script that you can download separately.
Each section includes an explanation followed by the corresponding
code.
Click the Code button in the bottom right corner to
display the code for each step/plot we are making.
Before we start, we need to install and load the required packages to
run the script.
The code below creates a list of the names of all packages we
need.
The next part checks if you have these packages installed, and if not,
installs them—before they are loaded into R.
# List of required packages
packages <- c(
"ggplot2","dplyr","tidyr","forcats","scales",
"viridis","RColorBrewer","ggsci","ggridges","ggdist"
)
# Install any packages that are not already installed
new.packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new.packages)) install.packages(new.packages)
# Load all the packages without printing the output
invisible(lapply(packages, library, character.only = TRUE))The default theme of ggplot has a gray background designed to make graphical elements “pop” out. For a cleaner look, we define our own custom theme that we will use for all figures. By defining theme_nice(), we can use this function in all our plots by adding + theme_nice(), saving time and ensuring consistency.
To easily add dimension and response level labels/names to our plots, we create vectors of dimension labels and response level labels. Below we specify both English and Norwegian versions, so you can switch between them as needed.
# Dimension labels (English & Norwegian)
eq5d_dim_eng <- c(
MO = "Mobility (MO)",
SC = "Self-Care (SC)",
UA = "Usual Activities (UA)",
PD = "Pain/Discomfort (PD)",
AD = "Anxiety/Depression (AD)"
)
eq5d_dim_nor <- c(
MO = "Mobilitet (MO)",
SC = "Egenomsorg (SC)",
UA = "Vanlige aktiviteter (UA)",
PD = "Smerte/ubehag (PD)",
AD = "Angst/Depresjon (AD)"
)
# Response level labels
eq5d_levels <- c("1","2","3","4","5")
eq5d_labels_eng <- c(
"No problem","Slight problem","Moderate problems",
"Severe problems","Extreme problems"
)
eq5d_labels_nor <- c(
"Ingen problemer","Litt problemer","Moderat problemer",
"Alvorlige problemer","Ekstreme problemer"
)If you have your own EQ-5D dataset, you can use the code below to load it here. Your data should be in standard wide format—each EQ-5D dimension must be stored as its own variable (column). Each dimension should be numeric, with levels ranging from 1 = No problems to 5 = Extreme problems. The EQ_VAS variable should be numeric, with a theoretical range from 0–100.
| Variable Name | Full Domain Name |
|---|---|
MO |
Mobility |
SC |
Self-Care |
UA |
Usual Activities |
PD |
Pain/Discomfort |
AD |
Anxiety/Depression |
EQ_VAS |
Visual Analogue Scale (0-100) |
Below are commented code examples for loading data from different file formats. Uncomment the relevant lines and modify the file path as needed.
# For SPSS (.sav):
# library(haven)
# df <- read_sav("your_data.sav") %>%
# mutate(across(where(is.labelled), as_factor))
# For Stata (.dta):
# library(haven)
# df <- read_dta("your_data.dta")
# For Excel (.xls/.xlsx):
# library(readxl)
# df <- read_excel("your_data.xlsx")
# For CSV (.csv):
# library(readr)
# df <- read_csv("your_data.csv")In order to have some data to showcase the various ways to visualize the EQ-5D instrument, we simulate a dataset of 2,000 responses. This includes demographic variables, hospital effects, and a sine-wave year effect. If you’re using your own data, skip this section.
set.seed(123)
n <- 2000
# Demographics & groups
gender <- sample(c("male","female"), n, replace = TRUE)
age_numeric <- sample(18:100, n, replace = TRUE)
age_breaks <- c(18,30,40,50,60,70,80,90,101)
age_labels <- c("18-29","30-39","40-49","50-59","60-69","70-79","80-89","90+")
age <- cut(age_numeric, breaks = age_breaks, labels = age_labels, right = FALSE)
treatment <- factor(sample(c("pre","post"), n, replace = TRUE),
levels = c("pre","post"), ordered = TRUE)
marital_status <- factor(sample(c("married","divorced/separated","single","widowed"),
n, replace = TRUE, prob = c(0.5,0.2,0.2,0.1)),
levels = c("married","divorced/separated","single","widowed"),
ordered = TRUE)
smoking <- sample(c("yes","no"), n, replace = TRUE, prob = c(0.3,0.7))
education <- factor(sample(c("Low","Medium","High"), n, replace = TRUE, prob = c(0.3,0.4,0.3)),
levels = c("Low","Medium","High"), ordered = TRUE)
income <- factor(sample(c("Low","Medium","High"), n, replace = TRUE, prob = c(0.3,0.4,0.3)),
levels = c("Low","Medium","High"), ordered = TRUE)
# Hospital effect
hospital_levels <- paste0("hospital ", 1:15)
hospital <- factor(sample(hospital_levels, n, replace = TRUE),
levels = hospital_levels, ordered = TRUE)
hospital_offsets <- seq(-25,25,length.out=15) + rnorm(15,0,3)
names(hospital_offsets) <- hospital_levels
# Simulate EQ-5D dimensions (1–5)
simulate_dim <- function(mean_target, smoker, age, education) {
b <- round(rnorm(1, mean = mean_target, sd = 1))
b <- b + ifelse(smoker=="yes", 2, 0)
b <- b + ifelse(age %in% c("60-69","70-79","80-89","90+"), 1,
ifelse(age %in% c("18-29","30-39"), -1, 0))
b <- b + ifelse(education=="Low", 1, ifelse(education=="High", -1, 0))
round(min(max(b, 1), 5))
}
MO <- mapply(simulate_dim, 2.5, smoking, age, education)
SC <- mapply(simulate_dim, 2.8, smoking, age, education)
UA <- mapply(simulate_dim, 3.1, smoking, age, education)
PD <- mapply(simulate_dim, 3.4, smoking, age, education)
AD <- mapply(simulate_dim, 3.8, smoking, age, education)
# Simulate EQ VAS (0–100)
simulate_vas <- function(baseline, gender, income, education,
marital_status, treatment, smoker, age, hospital) {
b <- qnorm(runif(1, pnorm(-1.2), pnorm(1.2))) * 8 + baseline
b <- b + ifelse(gender=="male", 3, -3)
b <- pmin(pmax(b, 0), 100)
b <- b + ifelse(income=="Low",-15, ifelse(income=="High",15,0))
b <- b + ifelse(education=="Low",-10, ifelse(education=="High",7,0))
b <- b + ifelse(marital_status=="married",8,
ifelse(marital_status=="divorced/separated",-7,
ifelse(marital_status=="widowed",-10,0)))
b <- b + ifelse(treatment=="post",10,0)
b <- b + ifelse(smoker=="yes",-10,0)
b <- b - ifelse(age %in% c("60-69","70-79","80-89","90+"), 8, 0)
b <- b + hospital_offsets[as.character(hospital)]
round(min(max(b, 0),100))
}
EQ_VAS <- mapply(simulate_vas, baseline = 60,
gender, income, education, marital_status,
treatment, smoking, age, hospital)
# Assemble and add year + sine-wave effect
df <- data.frame(
MO, SC, UA, PD, AD,
gender, age, EQ_VAS, treatment,
marital_status, smoking, education, income,
hospital
) %>%
mutate(
year = rep(2018:2025, each = n()/8),
EQ_VAS = round(pmin(100, pmax(0,
EQ_VAS + sin(2*pi*(year-2018)/7)*10)))
)Here, we convert our dataset from a wide format (where each EQ-5D dimension is in its own column) to a long format. This transformation makes it easier to create faceted plots in ggplot2.
Below you will see how the data looks before and after pivoting:
df_long <- df %>%
pivot_longer(
cols = c(MO,SC,UA,PD,AD),
names_to = "Dimension",
values_to = "Response"
)
# View data before pivot
head(df)## MO SC UA PD AD gender age EQ_VAS treatment marital_status smoking
## 1 5 5 5 5 5 male 60-69 56 pre married yes
## 2 5 5 5 5 5 male 60-69 42 pre divorced/separated yes
## 3 3 4 4 5 4 male 50-59 46 pre single yes
## 4 3 4 4 5 3 female 90+ 53 pre married no
## 5 2 1 1 5 2 male 18-29 52 post single no
## 6 5 5 5 5 5 female 80-89 25 pre married no
## education income hospital year
## 1 Low Medium hospital 11 2018
## 2 Low High hospital 5 2018
## 3 High Low hospital 10 2018
## 4 Medium High hospital 3 2018
## 5 High Low hospital 2 2018
## 6 Low Medium hospital 3 2018
## # A tibble: 6 × 12
## gender age EQ_VAS treatment marital_status smoking education income hospital
## <chr> <fct> <dbl> <ord> <ord> <chr> <ord> <ord> <ord>
## 1 male 60-69 56 pre married yes Low Medium hospita…
## 2 male 60-69 56 pre married yes Low Medium hospita…
## 3 male 60-69 56 pre married yes Low Medium hospita…
## 4 male 60-69 56 pre married yes Low Medium hospita…
## 5 male 60-69 56 pre married yes Low Medium hospita…
## 6 male 60-69 42 pre divorced/sepa… yes Low High hospita…
## # ℹ 3 more variables: year <int>, Dimension <chr>, Response <dbl>
Having followed the steps above, we are now ready to create data visualizations. We will start by plotting the EQ-5D dimensions to investigate distributions of responses across:
Most of the plots we will create are facet plots, meaning each figure is subdivided into panels—one per dimension—so we can easily compare response patterns.
This section presents bar charts that visualize the proportion of respondents selecting each category.
How this works:
eq5d_labels_eng, and the facets
are labeled via eq5d_dim_eng.eq5d_labels_nor and
eq5d_dim_nor.group_by(Dimension, Response) and
summarise(count = n()) tally the number of observations in
each response category for each dimension.group_by(Dimension) +
mutate(proportion = count / sum(count)) converts counts
into proportions, so that within each facet the bars sum to 1 (i.e.,
100%).scale_x_discrete() is used to split the x variable
labels, to avoid overplottingdf_bar <- df_long %>%
group_by(Dimension, Response) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(Dimension) %>%
mutate(proportion = count / sum(count)) %>%
ungroup()
p_bar_all <- ggplot(df_bar,
aes(
x = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng),
y = proportion,
fill = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng)
)) +
geom_bar(stat = "identity") +
geom_text(aes(label = paste0(round(proportion*100,1), "%")),
vjust = -0.3, size = 3) +
facet_wrap(~ Dimension, labeller = labeller(Dimension = eq5d_dim_eng)) +
scale_y_continuous(labels = scales::percent_format(), expand = expansion(mult = c(0,0.05))) +
scale_fill_brewer(palette = "Blues", name = "", labels = eq5d_labels_eng) +
labs(title = "Distribution of EQ-5D responses", x = "", y = "Percent") +
scale_x_discrete(
labels = function(x) stringr::str_wrap(x, width = 8)
) +
theme_nice() +
theme(
axis.text.x = element_text(size = 10)
)
p_bar_allTo compare distributions across gender, we add gender to the grouping and map it to the fill aesthetic.
How this works:
group_by()ensures that proportions
are calculated separately for males and females within each
dimension.position = position_dodge() places bars side-by-side
instead of stacking.fill = gender automatically assigns each
gender a distinct color from the Chicago palette.df_bar_gender <- df_long %>%
group_by(Dimension, Response, gender) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(Dimension, gender) %>%
mutate(proportion = count / sum(count)) %>%
ungroup()
p_bar_gender <- ggplot(df_bar_gender,
aes(
x = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng),
y = proportion,
fill = gender
)) +
geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
geom_text(aes(label = paste0(round(proportion*100,1), "%")),
position = position_dodge(width = 0.9), vjust = -0.3, size = 3) +
facet_wrap(~ Dimension, labeller = labeller(Dimension = eq5d_dim_eng)) +
scale_fill_uchicago() +
labs(title = "EQ-5D responses by gender", x = "", y = "Percent", fill = "Gender") +
scale_x_discrete(
labels = function(x) stringr::str_wrap(x, width = 8)
) +
scale_y_continuous(labels = scales::percent_format(), expand = expansion(mult = c(0,0.05))) +
theme_nice() +
theme(
axis.text.x = element_text(size = 10)
)
p_bar_genderHow this works:
df_bar, which already contains the
precomputed proportion for each (Dimension, Response) pair.filter(Dimension == "MO") isolates only the Mobility
rows. Adapt to choose the dimension of interestaes(), x = factor(Response…) ensures
the five response levels are ordered and labeled
correctly,y = proportion uses our precomputed proportions,
and fill colors each bar by its response level.geom_bar(stat = "identity") draws bars with heights
equal to our proportion values (instead of letting ggplot count the
data).geom_text()adds the percentage labels just above each
bar.df_bar_mo <- df_bar %>% filter(Dimension == "MO")
p_bar_mo <- ggplot(df_bar_mo,
aes(
x = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng),
y = proportion,
fill = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng)
)) +
geom_bar(stat = "identity") +
geom_text(aes(label = paste0(round(proportion*100,1), "%")),
vjust = -0.3, size = 3) +
scale_fill_brewer(palette = "Blues", name = "", labels = eq5d_labels_eng) +
labs(title = "Mobility (MO) response distribution", x = "", y = "Percent") +
scale_y_continuous(labels = scales::percent_format(), expand = expansion(mult = c(0,0.05))) +
theme_nice()
p_bar_moHow this works:
df_long, filter to
Dimension == "MO", then
group_by(Response, gender) to get counts per response per
gender.proportion = count / sum(count) gives us the share of each
response within each gender.fill = gender assigns a distinct color to
each gender, and position = position_dodge() places their
bars side-by-side.geom_text() is also dodged so the percentage labels
align correctly above each gender’s bardf_bar_mo_g <- df_long %>%
filter(Dimension == "MO") %>%
group_by(Response, gender) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(gender) %>%
mutate(proportion = count / sum(count)) %>%
ungroup()
p_bar_mo_g <- ggplot(df_bar_mo_g,
aes(
x = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng),
y = proportion,
fill = gender
)) +
geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
geom_text(aes(label = paste0(round(proportion*100,1), "%")),
position = position_dodge(width = 0.9), vjust = -0.3, size = 3) +
scale_fill_uchicago() +
labs(title = "MO response distribution by gender", x = "", y = "Percent", fill = "Gender") +
scale_y_continuous(labels = scales::percent_format(), expand = expansion(mult = c(0,0.05))) +
theme_nice()
p_bar_mo_gStacked bar charts are a handy way to visualize the distribution of responses across all five dimensions, and to compare groups at the same time. They show the composition of each dimension’s responses as 100% stacks, and can be faceted and/or grouped.
How this works:
Response to a factor with proper labels so
the stacks appear in the correct order.group_by(Dimension, Response)(and +
smoking for the grouped version) aggregates counts.mutate(pct = count / sum(count) * 100) converts counts
into percentages for each Dimension (or Dimension + smoking).position = "stack", geom_bar() layers
the categories on top of each other to reach 100%.facet_wrap(~ Dimension) splits the plot into panels by
dimension.# Prepare data for the stacked bar chart
df_stack <- df_long %>%
mutate(Response = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng)) %>%
group_by(Dimension, Response) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(Dimension) %>%
mutate(pct = count / sum(count) * 100) %>%
ungroup()
# Create the stacked bar chart
p_stack_all <- ggplot(df_stack, aes(x = Dimension, y = pct, fill = Response)) +
geom_bar(stat = "identity", position = "stack", colour = "white") +
geom_text(aes(label = paste0(round(pct,1), "%")),
position = position_stack(vjust = 0.5), size = 3) +
scale_y_continuous(labels = function(x) paste0(x, "%"),
expand = expansion(mult = c(0,0.02))) +
scale_fill_brewer(palette = "Blues", name = "", labels = eq5d_labels_eng) +
labs(title = "Stacked EQ-5D Responses (no group)",
x = "Dimension",
y = "Percent") +
theme_nice()
p_stack_allHere we break down the stacked bars by smoking status, faceted by dimension.
Only change from previous plot:
smoking to the group_by()
argumentssmoking to the x = argument# Prepare data for the grouped stacked bar chart
df_stack_smoke <- df_long %>%
mutate(Response = factor(Response, levels = eq5d_levels, labels = eq5d_labels_eng)) %>%
group_by(Dimension, smoking, Response) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(Dimension, smoking) %>%
mutate(pct = count / sum(count) * 100) %>%
ungroup()
# Create the facet-wrapped stacked bar chart
p_stack_smoke <- ggplot(df_stack_smoke, aes(x = smoking, y = pct, fill = Response)) +
geom_bar(stat = "identity", position = "stack", colour = "white") +
#uncomment to add percentages as text
#geom_text(aes(label = paste0(round(pct,1), "%")),
# position = position_stack(vjust = 0.5), size = 3) +
facet_wrap(~ Dimension, labeller = labeller(Dimension = eq5d_dim_eng)) +
scale_y_continuous(labels = function(x) paste0(x, "%"),
expand = expansion(mult = c(0,0.02))) +
scale_fill_brewer(palette = "Blues", name = "", labels = eq5d_labels_eng) +
labs(title = "Stacked EQ-5D Responses by Smoking",
x = "Smoking Status",
y = "Percent") +
theme_nice()
p_stack_smokeThe EQ VAS score allows for creating a variety of plots, depending on interest. Here, the script will showcase how to look at distirbutions across different groups (histograms, densityplots, boxplots, ridgline plot), and how to examine linear (and curvelinear) relationships between variables of interest and EQ VAS
This histogram shows the overall distribution of EQ VAS scores.
How this works:
aes(x = EQ_VAS) maps the VAS scores to the x-axis.geom_histogram(binwidth = 5)groups data into bins of
width 5; adjust binwidth to refine granularity.fill = "grey" sets the colour of the bars, whereas
colour = "white colours the gaps between. Set the colour
you like.p_vas_hist <- ggplot(df, aes(x = EQ_VAS)) +
geom_histogram(binwidth = 5, fill = "grey", colour = "white") +
labs(title = "EQ VAS Histogram",
x = "EQ VAS Score",
y = "Count") +
theme_nice()
p_vas_histCompare the VAS distribution across education levels with side-by-side bars.
How this works:
fill = education assigns each level a color.position = "dodge" arranges bars side-by-side.scale_fill_uchicago().p_vas_hist_ed <- ggplot(df, aes(x = EQ_VAS, fill = education)) +
geom_histogram(position = "dodge", binwidth = 5, colour = "white") +
scale_fill_uchicago() +
labs(title = "EQ VAS Histogram by Education",
x = "EQ VAS",
y = "Count",
fill = "Education") +
theme_nice()
p_vas_hist_edOverlay density curves for each education level, with group-specific mean lines. How this works:
geom_density() sets the plot to a density plot, and
requires only an x variable specified.fill and colour to
education draws separate density curves by levels of
education.alpha) helps visualize overlaps.geom_vline() highlight group
meansp_vas_den_ed <- ggplot(df, aes(x = EQ_VAS, fill = education, colour = education)) +
geom_density(alpha = 0.5) +
scale_fill_uchicago() +
scale_colour_uchicago() +
labs(title = "EQ VAS Density by Education",
x = "EQ VAS",
y = "Density") +
theme_nice()
# Add group-specific mean lines
group_vas_means <- df %>%
group_by(education) %>%
summarise(mean_vas = mean(EQ_VAS, na.rm = TRUE), .groups = "drop")
p_vas_den_ed <- p_vas_den_ed +
geom_vline(data = group_vas_means,
aes(xintercept = mean_vas, colour = education),
linetype = "dashed", linewidth = 0.5)
p_vas_den_edBoxplots are especially handy when you need to compare many groups at once—they show medians, IQRs, and outliers compactly.
How this works:
aes(x = EQ_VAS, y = hospital) rotates the usual
orientation so hospitals run vertically.coord_cartesian(xlim = c(0,100)) forces the x-axis to
span the entire theoretical range of the EQVAS score.p_vas_box <- ggplot(df, aes(x = EQ_VAS, y = hospital)) +
geom_boxplot(colour = "black", fill = "#6495ED") +
coord_cartesian(xlim = c(0,100)) +
labs(title = "EQ VAS Boxplot by Hospital",
x = "EQ VAS Score",
y = "Hospital") +
theme_nice()
p_vas_boxRidgeline plots overlay multiple density estimates, making it easy to compare distribution shapes across many groups.
How this works:
aes(y = hospital) creates a separate density for each
hospital.fill = after_stat(x) applies a gradient based on the
x-value (VAS score).rel_min_height and scale control overlap
and vertical spacing.scale_fill_viridis_c(option = "viridis") creates a nice
colour gradient (enter ?scale_fill_viridis in console to
learn more)p_vas_ridge <- ggplot(df, aes(x = EQ_VAS, y = hospital, fill = after_stat(x))) +
geom_density_ridges_gradient(scale = 1, rel_min_height = 0.01) +
scale_fill_viridis_c(option = "viridis", name = "EQ VAS") +
labs(title = "EQ VAS Ridgeline by Hospital",
x = "EQ VAS Score",
y = "Hospital") +
theme_nice()
p_vas_ridge## Picking joint bandwidth of 6.11
We now examine trends over time using both linear and smooth fits.
Observed Mean EQ VAS by Year & Gender
Plot the annual mean VAS for each gender with lines and points.
How this works:
stat_summary(fun = mean, geom = "line") connects yearly
means.stat_summary(fun = mean, geom = "point") overlays mean
points.Color and fill mapping to gender
differentiates the series.p_trend_obs <- ggplot(df, aes(x = year, y = EQ_VAS, colour = gender, fill = gender)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", shape = 21, size = 5) +
scale_fill_uchicago() +
scale_colour_uchicago() +
labs(title = "Observed Mean EQ VAS by Year & Gender",
x = "Year",
y = "Mean EQ VAS") +
theme_nice()
p_trend_obsAdd linear regression lines and confidence bands; you can swap or customize the smoother.
How this works & Tweaks:
method = "lm"fits OLS regression lines; use method =
“loess” for local smoothing.formula = y ~ poly(x, 2)or
poly(x,3) etc. inside geom_smooth().se = FALSE to hide confidence intervals.p_trend_lm <- ggplot(df, aes(x = year, y = EQ_VAS, colour = gender, fill = gender)) +
geom_smooth(method = "lm", se = TRUE) +
scale_fill_uchicago() +
scale_colour_uchicago() +
labs(title = "Linear Trend EQ VAS by Gender",
x = "Year",
y = "EQ VAS") +
theme_nice()
p_trend_lm## `geom_smooth()` using formula = 'y ~ x'
Capture non-linear patterns over time with LOESS smoothing; adjust span to control wiggliness.
How this works & tweaks:
p_trend_loess <- ggplot(df, aes(x = year, y = EQ_VAS, colour = gender, fill = gender)) +
geom_smooth(method = "loess", se = TRUE, span = 0.75) +
scale_fill_uchicago() +
scale_colour_uchicago() +
labs(title = "LOESS Trend EQ VAS by Gender",
x = "Year",
y = "EQ VAS") +
coord_cartesian(ylim = c(0,100)) +
theme_nice()
p_trend_loess## `geom_smooth()` using formula = 'y ~ x'
Facet the LOESS trend by age group and color by education level. How this works:
facet_wrap(~ age) splits panels by age group.Color and fill mapped to education
differentiate trend lines within each panep_trend_loess_age <- ggplot(df, aes(x = year, y = EQ_VAS, colour = education, fill = education)) +
geom_smooth(method = "loess", se = TRUE, alpha = 0.2) +
facet_wrap(~ age) +
coord_cartesian(ylim = c(0,100)) +
scale_fill_uchicago() +
scale_colour_uchicago() +
labs(title = "LOESS Trend of EQ VAS by Education & Age",
x = "Year",
y = "EQ VAS") +
theme_nice()
p_trend_loess_age## `geom_smooth()` using formula = 'y ~ x'
(aes()): map variables to
x, y, fill, and
colour.group_by(), summarise(),
and mutate() to compute counts, proportions, or means
before plotting.stat="identity" ↔︎
stat="count", adjust position ("dodge",
"stack", "fill"), or recode factors for
order.binwidth, breaks, or
boundary to refine. Use position="fill" for
relative frequencies.bw, set
alpha for transparency.geom_smooth(), change method (e.g.,
"lm", "loess", "gam"), tweak
span, or supply formula for polynomial fits.facet_wrap(~ variable) or
facet_grid(rows ~ cols) for multi-panel layouts.scale_fill_brewer(), scale_fill_viridis_d(),
or scale_fill_uchicago(), or define custom palettes with
scale_fill_manual().ggsave() calls to export figures;
adjust width, height, and dpi to match your publication
requirements.